Astronomy <fc Astrophysics manuscript no. 0417arti 


2nd February 2008 


(DOI: will be inserted by hand later) 





Simulations of dust-trapping vortices in protoplanetary discs 

Anders Johansen 1 ' 2 , Anja C.Andersen 1 , and Axel Brandenburg 1 



^1" 

o 
o 

(N 
C 

a 
>—> 

ON 
(N 

> 

in 
o 
o 

m 
o 

^ ' 

6 ■ 
a: 

C/3 . 



1 NORDITA, Blegdamsvej 17, DK-2100 Copenhagen, Denmark 
e-mail: aj ohan/anj a/brandenb@nordita . dk 

2 Astronomical Obs., NBIfAFG, Copenhagen University, Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark 



2nd February 2008, Revision: 1.2 

Abstract Local three-dimensional shearing box simulations of the compressible coupled dust-gas equations are 
used in the fluid approximation to study the evolution of different initial vortex configurations in a protoplanetary 
disc and their dust-trapping capabilities. The initial conditions for the gas are derived from an analytic solution to 
the compressible Euler equation and the continuity equation. The solution is valid if there is a vacuum outside the 
vortex. In the simulations the vortex is either embedded in a hot corona, or it is extended in a cylindrical fashion 
in the vertical direction. Both configurations are found to survive for at least one orbit and lead to accumulation 
of dust inside the vortex. This confirms earlier findings that dust accumulates in anticyclonic vortices, indicating 
that this is a viable mechanism for planetesimal formation. 
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1. Introduction 

The expected scenario for planet formation within a proto- 
planetary disc around a newly formed star is that planets 
grow from kilometre-sized planetesimals to protoplanets 
by gravitationally attracting each other. The planetesi- 
mals are expected to form within the protoplanetary disc 
while gas is still present in the disc. The coupling of gas 
and solids is therefore an important issue in the formation 
of planets. 

Planetesimals can only grow from sub-micron solids if 
their relative velocities are less than about 1 ms _1 , since 
larger velocities will result in disruption of the aggregate 
(Blum & Wurm 2000). Relative velocities considerably 
higher than 1 ms _1 arise from chaotic motions in a turbu- 
lent disc or from differential orbital drift in a laminar disc 
(Weidenschilling & Cuzzi m?9"3|l . Another problem facing 
planetesimal growth is the rapid inward orbital drift asso- 
ciated with gas drag that carries small bodies into the star 
on a time scale of 100-1000 years (Weidenschilling Il977fl . 

One attractive scenario which may be able to over- 
come the two above problems is the presence of dust- 
trapping anticyclonic vortices within the protoplanetary 
disc (Barge & Sommeria 1995; Tanga et al. 119961 Bracco 
et al. RUM Godon & Livio |2TjOTJ|> . If dust is trapped 
within vortices, the relative velocities between solid parti- 
cles would be small since all the solids rotate in the same 
direction within the vortex. Trapping would also prevent 



solid bodies from spiralling inwards towards the star. Self- 
gravity between sufficiently large amounts of trapped solid 
material inside vortices could even trigger a local gravi- 
tational instability and subsequent growth of centimetre- 
sized solid bodies into planetesimals. 

The formation and stability of vortices has been an 
area of much research since the dust-trapping mechanism 
was proposed by Barge & Sommeria 11995} . Using a two- 
dimensional incompressible model, Bracco et al. ([1998 
1999} show that anticyclonic vortices can form as a relic 
from the originally strongly turbulent disc. The existence 
of a baroclinic instability in a protoplanetary disc is inves- 
tigated by Klahr & Bodenheimer |2003). This instability 
works as a source of vorticity and forms vortices similar 
to how vortices form around high and low pressures on 
Earth. The stability of vortices is simulated by Godon 
& Livio (199SJ. They show that two-dimensional anticy- 
clonic vortices can survive for hundreds of orbits, and that 
isolated vortices can merge to form larger vortices. 

The dust-trapping mechanism was investigated by 
Hodgson & Brandenburg l|1998[) in three dimensional 
simulations of disc turbulence driven by the magneto- 
rotational instability (Balbus & Hawlev Il998f) . In these 
simulations the vortices are highly time dependent, but 
when gas velocity is frozen in time, they find concentra- 
tion of dust inside vortices. However, when considering a 
freely evolving gas flow, they no longer see any significant 
dust concentration inside vortices. This is associated both 
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with the finite life time of turbulent vortices and with the 
stirring up of dust by turbulence in the vertical direction. 
The dust-trapping efficiency of vortices was explored an- 
alytically by Chavanis ( 2000 ) for vortices of arbitrary as- 
pect ratio. Recently de la Fuente Marcos & Barge ( 2001) 
have considered a single frozen vortex velocity field in two 
dimensions and a distribution of particle sizes. Using re- 
alistic expressions for the drag force on dust they find 
efficient particle trapping inside vortices, preventing the 
inwards orbital drift associated with gas drag. 

In the present paper we explore the dust-trapping 
efficiency of three-dimensional vortices. We start out 
from the analytic vortex solution by Goodman et al. 
( |1987, hereafter referred to as GNG) ). This solution takes 
density and pressure gradients across the vortex explicitly 
into account. We let the gas flow evolve freely and exa- 
mine the coherence of the three-dimensional vortices and 
their ability to trap dust. 

The paper is structured as follows. In Sect. we give 
the basic hydrodynamic equations used. In Sect.Olwe dis- 
cuss in detail the vortex solution that we use as initial 
condition. The dust-trapping mechanism is described in 
Sect.0] In Sect.[5]we discuss the two models that we use. In 
Sect. we describe the numerical scheme and the bound- 
ary conditions implemented. Our results are presented in 
Sect.EI 

2. Basic hydrodynamic equations 

We perform simulations in the shearing sheet approxi- 
mation (Wisdom & Tremaine 119881 Hawley et al. I1995JI . 
where a local coordinate frame corotating with the Kepler 
flow at a distance ro from the central star is considered. 
In this local approximation the curvilinearity of the coor- 
dinates is neglected, so the validity is limited to the case 
when the size of the vortex is much smaller than ro. The 
x-axis points away from the star, and the y-axis points 
along the flow direction. The angular velocity profile of a 
Keplerian disc goes as fi(r) oc r~ 9 , where q = 1.5 when 
considering only the gravitational attraction to the central 
star. 

In the coordinate frame both inertial and fictitious 
forces exist. Radial gravity and centrifugal force terms 
only cancel at the radius ro, giving rise to a tidal force 
away from ro- However, when measuring velocities rela- 
tive to the main shear flow, = (0, —qQx, 0), this tidal 
force vanishes. We then have u = + u and use u as 
the velocity variable. The vertical part of the gravity still 
gives rise to a restoring force proportional to — f2 2 z in the 
z-direction. Gas particles also experience a pressure gra- 
dient force, viscous forces and the fictitious Coriolis force. 

We treat gas and dust as two separate fluids and let 
the two interact through a drag force. We use the term 
dust for solid bodies of all sizes. The drag force from gas 
upon the dust attempts to accelerate dust to match the 
velocity of the gas, and vice versa. We assume that dust 
is not pressure supported (due to a low number density 
and low particle velocities). Since the solid density of dust 



particles is many orders of magnitude higher than the gas 
density and the pressure gradient acts as a volume force, 
it is reasonable to assume that dust is not affected by gas 
pressure (see e.g. Seinfeld [$86). 

2.1. The shearing sheet approximation 

In the shearing sheet approximation the equations for the 
departure from the main shear flow (Brandenburg et al. 
1995) can be written in the form 
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pTl^+u- Vsj = 2pvS 2 +pQ{V-u) 2 + V-{KVT) , (5) 

where u, P, p, T and s are the velocity, pressure, density, 
temperature and specific entropy of gas, respectively, and 
u d and pd are the velocity and density of dust. The density 
of dust is defined as p d = n d m d , where n d is the number 
density of dust particles and m d is their mass (all particles 
are assumed to have the same mass). The parameter £ de- 
scribes some bulk viscosity, see Sect. 12.21 Dust density is 
not to be confused with the solid density p s of individual 
dust particles (defined in Sect. 12. The advective deriva- 
tive, V /Vt = d/dt+Uyd/dy, is with respect to the mean 
shear flow only, 



IflUy 



(6) 



z 



is the Coriolis force combined with the tidal force and the 
vertical gravity force, F v i sc and -F v isc,d are viscosity forces 
(see next section) and (3 is the coupling coefficient between 
gas and dust. This coupling coefficient can be expressed 
in terms of the stopping time r s , which in turn is defined 
through the relation 



Pd 

(Md 



(7) 



where Fd is the drag force per unit volume, so /3 — p d /r s . 
The stopping time is a parameter that describes the inter- 
action between a single dust particle and the surrounding 
gas and does not depend on the number density of dust 
particles, whereas (3 does. 

Temperature, pressure and density are related to each 
other by the relation P = (c p — c v )Tp, where c p and c v 
are the specific heats at constant pressure and constant 
volume. The adiabatic sound speed c s is given by 



<4 ex P bs/cp + (7 - 1) ln(p/p )] , 



(8) 
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where po and c s o are arbitrarily chosen integration con- 
stants from the integration of the first law of thermo- 
dynamics and 7 = c p /c v is the ratio of specific heats at 
constant pressure and volume, respectively. We have cho- 
sen integration constants such that s = when c s = c s0 
and p = pq. The thermal conductivity K is related to the 
kinematic (shear) viscosity v through the non-dimensional 
Prandtl number, Pr = v/(pc v K). 

2.2. Viscosity 

The viscosity force on the gas can be expressed as 
—^-=v (V 2 « + jVV ■ m + 2S • V lnp) + (VV • it, (9) 



where v is the kinematic viscosity of the gas (assumed 
constant), 



1 fdi 



8Un 



2 \ dxj dxi 



(10) 



is the traceless rate-of-strain tensor and C is the bulk vis- 
cosity, which is invoked solely to smear out sharp gra- 
dients (or shocks) over a few mesh points near strongly 
converging regions. The shock viscosity, Q — Cshock, as it 
is used here, is proportional to the smoothed maximum of 
the positive flow convergence, 



Cshock = c s hock ^max[(- V • u)+] ) , 



(11) 



where c s hock is a non-dimensional coefficient defining the 
strength of the shock viscosity. The smoothing is accu- 
rate to second order, and the maximum is taken over the 
second-nearest points. 

The viscous force on the dust would normally be neg- 
ligibly small, but when treating dust as a fluid one always 
has to add a small viscosity to damp out unphysical wig- 
gles on the mesh scale. It proved sufficient to use 



Pd 



1 F 

-*- v 



(12) 



where we have ignored any dependence on dust density, 
because -F v isc d is small anyway. 



2.3. The stopping time 

Due to Newton's third law the value of (3 must be the 
same for dust and gas, which leads to 



Pd 



P 



(13) 



Specifying the stopping time of dust automatically sets 
the stopping time of gas, r Sig , according to Eq. I|13fl. The 
stopping time of gas decreases with increasing dust den- 
sity. This means that the back-reaction of dust upon gas 
should become more and more important as dust density 
approaches gas density inside a vortex. 

A non-dimensional measure of the stopping time in 
terms of the angular velocity is the parameter i?r s . If (1t s 



is small, then dust is strongly coupled to gas, whereas 
when Qt s is large, dust is relatively unaffected by gas drag. 

Analytic expressions for the stopping time of dust due 
to gas drag exist in various regimes, depending on the 
mean free path of the particle and the Reynolds number of 
the flow (e.g. Weidenschilling 19771 see also Chavanis 2000 
de la Fuente Marcos & Barge ;2001j. The dust particles 
are assumed spherical with radius a s , solid density p s and 
mean free path A. If A < |a s , the Stokes drag law is valid. 
Here the stopping time depends on Reynolds number Re. 
When A > |a s , one must use the Epstein drag law. The 
linear drag law used here, where r s is constant in Eq. Q, 
is only valid when the stopping time is independent of 
relative velocity. This is the case in two regions: in the 
Stokes drag regime with Re < 1 and in the Epstein drag 
regime (see Weidenschilling H977fl . The transition between 
Epstein and Stokes regime in the Solar nebula depends on 
dust particle radius (Chavanis 2000). In the outer Solar 
System particles will typically be in the Epstein regime, so 
the drag law used here is valid for the outer Solar System, 
i.e. from the location of Jupiter and outwards. We ignore 
the dependence of r s on local gas density. 

In the Epstein drag regime there is a simple expression 
for the stopping time, 



P Cg 



(14) 



To calculate the radius of a particle whose stopping time 
is known, we consider a typical protoplanetary disc with 
a scale height H = 10 12 cm and a mass ratio p s /p — 10 10 
in the mid-plane. The scale height of the disc can be ex- 
pressed in terms of the speed of sound and the Kepler fre- 
quency under the assumption of vertical hydrostatic equi- 
librium and an isothermal density profile, 



H = c s /n, 



so 

P 2 

a s = — H[2t s = 10" S7t s cm . 

Ps 



(15) 



(16) 



2.4. Settling of dust 

Dust is subjected to a gravity force in the vertical direction 
without a balancing pressure gradient force, which makes 
it fall towards the mid-plane. The terminal velocity (the 
velocity where gravity and drag force balance) at height z 
is determined by the stopping time as 



9 1 

= -Q 2 z u 



(17) 



so = ~ T s^ 2 z. At z = H the terminal speed reaches 

the fraction t s (2 of the sound speed, in which case the 
linear drag law Eq. J7J) breaks down for bodies with large 
stopping times. 

For long stopping times, t s > Q , dust will essen- 
tially free-fall towards the mid-plane. The solid particles 
then perform damped oscillations around the mid-plane 
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and eventually settle to form a thin sheet around z — 0. 
This was long considered the seed of planetesimals: the 
thin sheet will become gravitationally unstable and then 
planetesimals will condense out (Safronovfl969. Goldreich 
& Ward IT§73J) . However, today it is believed that tur- 
bulence caused by the shear between the thin dust layer 
and the gas leads to a continuous stirring up of the dust 
(Weidenschilling 119801 Cuzzi et al. ^993 ) . In our present 
model the first arrival of the solid particles at the mid- 
plane (forming essentially a delta-function) leads to the 
dust density being under-resolved for given resolution. We 
avoid this problem by ignoring the vertical gravity on the 
dust thereby maintaining the initial scale height. 



3. Initial conditions 

3.1. The vortex solution 

GNG showed that there exists an elliptic vortex flow solu- 
tion to the Euler and continuity equations in the shearing 
sheet approximation. In this section we briefly describe 
their solution. In Appendix A we go into more detail about 
how to arrive at this solution. 

For a better distinction of the different forces involved 
we look at the equations with the main shear flow included 
again. The Euler and continuity equations for the motion 
of gas then take the explicit form (when ignoring drag 
force and viscosity) 



— + (u-V)u = n 2 (3x 



at 



2f2xu--VP, (18) 
P 

(19) 



where x = (x, 0, 0) and z = (0, 0, z). The 3f2 2 x term is the 
tidal force approximated along the x-direction to first or- 
der, and the — fl 2 z term is the vertical gravity force. The 
last two terms on the right hand side of Eq. (|18|) are the 
Coriolis force and the pressure gradient force. We stress 
that the equations in this form are completely similar to 
the Euler and continuity equations of Sect. 12.11 We intro- 
duce the specific enthalpy h and write — p~ 1 'VP = — V/i 
in the isentropic case. Then the GNG solution for the el- 
liptical velocity field and the corresponding enthalpy is 



u x 

Uy 

u z 

h = ^n*(b 2 -x A -eV)- 



fi'x . 

e 

0, 
1 



Q 2 z 2 . 



(20) 
(21) 
(22) 
(23) 



where e = a/b is the aspect ratio of the ellipse and Q' is 
the angular velocity of the vortex. The ellipse has a semi- 
major axis a in the y-direction (along the Kepler flow) and 
a semi- minor axis b in the x-direction (radially outward). 
The aspect ratio must be in the interval ^ e ^ 0.5. The 
velocity field is incompressible and has no z component. 
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Figure 1. The angular velocity of a vortex ft' as a func- 
tion of aspect ratio e. The dash-dotted line connecting the 
end-points is present to illustrate the curvature. 

The parameters fi 1 and 8 are connected to the aspect ra- 
tio of the ellipse and the background Keplerian angular 
velocity through 



n> = a Q = ^ a 
(see Fig. P), and 



8 2 = - 



2\/3 



1 



(24) 



(25) 



These two bindings are necessary for fulfilling the conti- 
nuity equation (see Appendix A). The solution is valid 
where h ^ 0, giving the vortex an ellipsoidal shape with 
axes a — b/e, b and c = 5b. 

The anticyclonic vortex feels a compressive Coriolis 
force, which is balanced by pressure and tidal forces. The 
anisotropy of the latter gives rise to the elliptical shape. In 
the pressureless limit, the streamlines reduce to Keplerian 
epicycles. We explain the vortex dynamics and the dust- 
trapping mechanism in Sect.^J 

3.2. Units 

The model is scale- invariant. This means that the basic 
units can be chosen arbitrarily. We take 



[t] = n-\ [ x ] = b, [p] = Po , 



(26) 



where b is the semi-minor axis of the vortex in the hori- 
zontal plane, i.e. in the cross-stream direction, and po is 
the average density of gas in the whole box. The latter is 
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a conserved quantity, since there is no flow through the 
boundaries. The unit of velocity is derived from the basic 
units as [u] = [x]/[t] = bfl, and the unit of acceleration 
is [a] = [x]/[t] 2 = bfl 2 . For clarity we will often write the 
units out explicitly. 

3.3. Global solutions 

The velocity field given by Eqs. 1(20 )1 -lf22 ]l is not a global 
solution, since the velocity field is discontinuous when 
crossing the vortex boundary to the surrounding Kepler 
flow, where the only velocity component is u y — —qflx. 
By looking at the y-velocity of the flow, 



-afJx 



V3 



-.fix . 



(27) 



it is apparent that, regardless of the choice of e, the tan- 
gential part of the vortex flow will always be faster than 
the Kepler flow at any position in x. It seems that for 
the vortices considered here there is no way of producing 
a linear velocity field that can make a gradual transition 
from the vortex to the surrounding Kepler flow. Our ini- 
tial conditions do therefore not satisfy our equations at 
the interface between vortex and exterior. For this reason 
we must expect to see an evolution away from our initial 
state. 
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Figure 2. Forces affecting the gas. Average magnitudes of 
pressure gradient force F p , Coriolis force Fcor plus tidal 
force F g+C and resulting centre-directed force F cen as func- 
tions of aspect ratio e. Averages are taken over the vortex 
boundary. The force unit is bfl 2 . 



4. Dust-trapping mechanism 

4.1. Forces driving the vortex 

The fact that anticyclonic vortices suck in dust particles 
can be explained by looking at the forces involved in the 
rotation. The rotation of gas is maintained because the 
Coriolis force F'cor, the tidal force -Fg+c and the pressure 
gradient force F p add up to a resulting force F ccn directed 
towards the centre of the coordinate system. The symbol 
F is here used for the force per unit mass. 

The pressure gradient force is calculated from Eq. (|23|l 
(we consider z = since the pressure gradient in the in- 
direction is completely balanced by the vertical gravity 
component) , 



-Vh 



5 2 f2 2 x 
5 2 fl 2 ^ 2 



e y 



2V3 



s 2 f> 2 



e 2 y 



1 



VT 
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(28) 



and the Coriolis force is 
F C or = -2J?fc x u = -Ifl 




= -2Q 



af2x/e S \ , } 




e 2 y 



2V3 



:Q 2 



e 2 y 



(29) 



Note that both forces point perpendicular to the vortex 
flow. The resulting force, 



F 




n 2 I y I , (30) 



points always towards the centre and thus produces the ro- 
tating motion. Fig.|2]is a plot of the average values over the 
vortex boundary of the pressure gradient force, Coriolis 
force plus tidal force and resulting centre-directed force as 
a function of e. The force magnitudes have the unit of bfl 2 . 
As the aspect ratio decreases, the pressure gradient force 
becomes relatively more important. At e = 0.5 the vortex 
is in equilibrium without any need for a pressure gradients 
in the x-y plane. This pressure-less vortex is exactly the 
epicyclic vortex considered by Barge & Sommeria (1995). 
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Figure 3. Forces affecting dust as a function of / (dust speed as a fraction of gas speed). Average values of Coriolis 
force Fcor plus tidal force F g+C , necessary centre-directed force to maintain rotation F con and excess force F cxc for 
different values of aspect ratio e. The force unit is bfl 2 . Averages are taken over the vortex boundary. The forces do 
not have the same direction, so the two vectors have been added to get the magnitude of the excess force. The excess 
force converges to the missing pressure gradient force for / = 1 in all cases. 



4.2. Forces on the dust 

Dust, initially at rest, is accelerated by drag to follow 
gas around the vortex, but in the beginning with a much 
smaller velocity. Moving at a fraction / of the gas speed 
relative to the main shear flow, u& = fu, we find for 
Ud = u<i + it(°) the expression = fu + (1 — f)u^°\ so 
dust is subjected to the Coriolis force 

-^f/n- ( A) - (i - /)(*?*). PD 

The centre-directed force needed to keep dust spinning 
around the vortex is proportional to velocity squared, and 



is therefore reduced by a factor of f 2 compared to the 
force needed to keep the gas rotating, 

of2 2 fx \ 
F ccn = ~Y^n 2 \y\. (32) 

The excess force on the dust is then -F e xc = F'cor + 
.Fg+c — F ccn . In Fig. |3 Coriolis force plus tidal force, 
centre-directed force needed to maintain the same radius 
of rotation as the gas and excess force working on the 
dust are shown as a function of / and for different values 
of the aspect ratio e. For all value of e the Coriolis force 
grows much faster with / than the required centre-directed 
force. This leads to an excess force directed inwards, which 
causes dust to spiral inwards. Similar explanations for the 
dust-trapping mechanism of vortices are given in Tanga et 
al. (1996) and Chavanis (2000). Already at a few percent 
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of the gas speed, the excess force on the dust is signifi- 
cant. Cyclonic vortices have an outwards directed Coriolis 
force. They will therefore not be able to concentrate solid 
particles, but will rather expel them. 

The centre-directed force catches up with the Coriolis 
force at / = 1.0 for e = 0.5, but for higher e there is still an 
excess force on the dust at f — 1.0 of the same magnitude 
(but opposite direction) as the pressure gradient force af- 
fecting the gas. The reason for this is that dust particles 
are not subjected to the pressure gradient force of gas, so 
even if dust particles were light enough to quickly become 
accelerated to match gas velocity, they will feel an extra 
force inward and thus begin to spiral inwards. The ter- 
minal velocity that dust obtains perpendicular to the gas 
flow, , is a result of a balance between excess inwards 
force and drag force, 



-5 2 Q 2 



e 2 y 



— u 



d_L 



(33) 



which means that 



never become negative, so embedding the vortex in a disc 
of finite density and constant entropy is not possible in 
the vertical direction. 

We use two different ways of modelling the vortices in 
three dimensions: hot corona and cylindrical vortex. 

5.1. Hot corona 

Hydrostatic equilibrium in the z-direction can be obtained 
by embedding the vortex in a tenuous gas of high tem- 
perature, i.e. a hot corona. In terms of specific enthalpy, 
Eq. I|37|l can be written as 



d<P 

dz 



h ds 
c p dz 



dh 

dz 



(39) 



For a given entropy distribution it is then possible to 
construct a continuous enthalpy so that hydrostatic equi- 
librium is obtained. This method is also used by von 
Rekowski et al. jl2UU51l . 

The hot corona simulations are done in a box of size 
(L x ,L y ,L z )/b= (4,32,4). 



,s 2 n 2 



e 2 y 



(34) 5.2. Cylindrical vortex 



in analogy to the vertical settling of dust towards the mid- 
plane due to gravity. The time scale for this is quite long, 



t 



1 



M 2 



(35) 



considering that the stopping time must be short for the 
dust to match the velocity of the gas, but if vortices are 
indeed long-lived, it could prove an important mechanism 
for trapping dust particles of short stopping time inside 
vortices of low aspect ratio. 

The aspect ratio also has an influence on dust-trapping 
at low /. It is evident from Fig. [3] that the excess force 
grows faster with / and reaches a higher maximum for 
high e. 

5. Three-dimensional models 

The vortex solution of GNG has zero enthalpy outside the 
vortex, corresponding to zero density through the relation 



Po 



(7-1) -re 



L -s() 



1/(7-1) 



(36) 



This gives rise to a potential problem in three dimensions: 
hydrostatic equilibrium in the z-direction requires 



d$ 1 dP _ ( | 

dz p dz 



where 
1 

2' 



<p = -n 2 z 2 



(37) 



(38) 



is the vertical gravitational potential. This means that the 
pressure must fall off in the vertical direction, but it must 



Here we extend the vortex to the entire z-length of 
the box. The vortex is then no longer ellipsoidal, but 
rather cylindrical with an ellipse-shaped cross section. 
Hydrostatic equilibrium is obtained without use of entropy 
through 



dh 

dz 



d<P 

dz 



(40) 



so h{z) — — ii? z + ho, where ho must be greater 
than ^f2 2 z 2 Qp to avoid negative enthalpy anywhere. The 
cylindrical vortex simulations are done in a box of size 
(L x , L y , L z )/b — (4,32,0.4). The reason for using a shal- 
lower box for the cylindrical vortex is to allow for a smaller 
ho and thus to have a large density ratio between the vor- 
tex and its surroundings (see Sect. 17.31) . The minimum en- 
thalpy in the mid-plane is ^fi 2 z 2 op , so lowering z top and 
adopting the minimal enthalpy possible thus leads to a 
larger density ratio. 

In both cases we let dust start with zero initial velocity. 
Dust density at height z is initialised to be a fraction r = 
0.01 of gas density (see Natta et al. l2000|l far away from the 
vortex at the same height. Although the actual value of r 
does not directly influence dust dynamics, it does influence 
the amount of back-reaction from dust upon gas. 

6. Numerical method and boundary conditions 

We simulate the motion of gas and dust in a coordinate 
frame corotating at the local Keplerian speed. We use the 
Pencil-Code 1 (Brandenburg & Dobler I2002[) which uses 
third order Runge-Kutta time stepping and a sixth order 



1 The code is available at 
http: //www.nordita. dk/data/brandenb/pencil-code 
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Figure 4. The initial condition in the mid-plane for the 
hot corona vortex: gas density and gas velocity (left) and 
the corresponding dust density and velocity (right). The 
velocity field does not seem to follow the contours of the 
vortex, but this is because the Kepler velocity has been 
subtracted. 



finite-difference scheme in space and employs central finite 
differences, so the extra cost of recentering a large number 
of variables between staggered meshes each time step is 
avoided. The code solves the non-conservative form of the 
equations (see Brandenburg 2003 for details). 

Periodic boundary conditions for velocities, density 
and specific entropy are used in the x- and y-directions. 
The periodic boundary condition in x is appropriately 
sheared in the shearing sheet approximation. In the z- 
direction we use stress-free boundary conditions, i.e. 



Udx, 



Udz 



0. 



(41) 



where commas denote partial differentiation. Dust and gas 
densities on the boundaries are fully determined by the 
equations, but in practice we need to specify values in 
the ghost zones which is accomplished by extrapolation. 
We use a resolution of (n x ,n y ,n z ) = (128,256,128) grid 
points. 




□ 
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Figure 5. The initial condition in the vertical for the hot 
corona vortex: gas density (left) and the corresponding 
dust density (right). Dust density as a function of z is 
initially 0.01 times gas density in the same height. The 
density colour scale is the same as in Fig. 0] 



7. Results 

7.1. Initial conditions and choice of stopping times 

The initial condition for gas and dust in the hot corona 
model is plotted in Fig. 0] (mid-plane , only one fourth of 
the y-length of the box is shown) and Fig. (x-z plane). 
We choose an aspect ratio of e = 0.4 since we expect dust- 
trapping to be most efficient for vortices of high aspect 
ratio. The transition from vortex to surroundings has been 
smeared out both in density and velocity field. We have 
experienced that the simulations run better when abrupt 
transitions are smoothed out. The dust is initially at rest. 
In all plots of velocity fields we use the Kepler-subtracted 
velocities u and ttd- 

We have run hot corona simulations for three different 
values of the stopping time: J7r s = 0.1 (dust completely 
coupled to gas), Qt s = 1.0 (dust semi-coupled to gas) 
and J?t s = 10.0 (dust almost not coupled to gas). This 
corresponds to rocks of sizes respectively a s = 10 cm, 
a s = 100 cm and a s = 1000 cm, respectively. For the 
minimum-mass nebula of Cuzzi et al. ( 1993) the mean free 
path is A = (r/AU) 11 / 4 cm. This means that we must, 
strictly speaking, go beyond r = 10AU for the largest 
particles to be in the Epstein regime. For a s — 10 cm 
and a s = 100 cm the Epstein drag law is valid already 
before r — 5 AU. The range of dust radii considered is 
expected for the largest agglomorations of sticking dust 
particles after having reached the mid-plane, according to 
Weidenschilling & Cuzzi (1993) who model how dust par- 
ticles falling towards the mid-plane stick to each other in a 
turbulent protoplanetary disc. Larger sizes are prohibited 
by the turbulence in the disc. 

The cylindrical vortex was only run for Qt s — 1.0, 
since we are mostly interested in whether it is a valid 3D 
vortex solution. Wc follow the evolution of gas and dust 
for a full Kepler orbit i max — 2n/f2 (which is not quite a 
full rotation of the vortex, see Fig.^). 
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Figure 6. Gas and dust in the mid-plane for the hot corona after one orbit. Gas density and velocity field for the three 
values of Qt s are indistinguishable, so only J?t s = 1.0 is shown. The other three plots show dust density and dust 
velocity field after one orbit for Qt s = 0.1, J?r s = 1.0 and Qt s — 10.0 respectively. The dust velocity field of Qt s = 0.1 
is very similar to that of the gas, due to the short stopping time. For fir s = 1.0 there is a strong convergence towards 
the interior of the vortex, whereas for J7r s = 10.0 only a slight density increase in a narrow region that extends from 
the vortex along the shear is seen. The density colour scales are the same as in Fig. 0] 



7.2. Viscosity parameters 

For the hot corona model we are able to run for a whole 
orbit with a viscosity of v = i/<j = 2 • 10 -4 and a shock vis- 
cosity of c s hock = 2. This corresponds to a mesh Reynolds 
number of 



Thus, since Re mes h is rather large, the viscosity v is al- 
most everywhere unimportant, and most dissipations oc- 
curs through the shock viscosity (not included in the ex- 
pression for Re mcs h), but this affects only convergent flow 
regions and not the vortical parts of the flow. This can be 
seen from the following consideration. 



Re m csh = max(|ii|) max(<5:c, Sy, 5z)/v ~ 625 . 



(42) 
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Figure 7. Gas and dust in the x-z plane for the hot corona case after one orbit. Gas density and velocity field for 
the three values of fir s are indistinguishable, so only Dt s = 1.0 is shown. The other three plots show dust density 
and dust velocity field after one orbit for J2t s = 0.1, Qt$ = 1.0 and Qt s = 10.0 respectively. For Qt s = 0.1 dust is so 
strongly coupled to gas that the velocity fields are very similar, whereas i?r s = 1.0 shows strong convergence of dust 
in the interior of the vortex. For Qt s — 10.0 dust is not as affected by gas motion, although strong motions move dust 
around the vortex and upwards. Velocities are exaggerated about 10 times compared to Fig. The density colour 
scales are the same as in Fig. 



The Laplacian of u can be written V 2 it = —V x V x 
u + VV • u. Using this we can rewrite the expression for 
the viscous force as 

= v (- V x u + 2S ■ V lnp) + (%v + Q VV • it , (43) 

P 

where u> = V x u is the vorticity of the flow. This shows 
that the kinematic shear viscosity works as a drain for 
vorticity in regions where the vorticity changes, whereas 
shock viscosity only affects regions of convergence in the 
velocity field. 

At the boundary between vortex and surroundings the 
vorticity changes abruptly, and we have experienced vor- 
tices dissolving rapidly when using too high a viscosity. 
We stress that this is a purely numerical issue, and that 
real physical viscosity is many orders of magnitude lower 
than the viscosity we use here. 

For the cylindrical vortex we had to use a viscosity of 
v = v& = 10~ 2 and a shock viscosity of c s hock = 4. The 
reason for these higher viscosities is that the cylindrical 
vortex does not have as large a density ratio as the hot 
corona vortex, so it takes more viscosity to keep the vortex 
intact on the boundary. 

7.3. Lifetime of vortices 

When using a realistic disc background density we have 
experienced the vortices breaking up at the local sound 
speed at the vortex boundary. This we believe is caused by 
the term V-iZ in the continuity equation (J3J. Although the 
initial condition has V • u — analytically, this is not nec- 
essarily valid numerically. Calculating the spatial deriva- 
tives of the velocity field on a Cartesian coordinate grid 
yields V • u ^ over a few mesh-points (due to the order 
of the numerical derivatives) around the vortex boundary. 
The result is the depletion of density at the NW and SE 



corners of the vortex, and increase of density at the NE 
and SW corners. We find that we can alleviate this prob- 
lem by having a large density contrast between the vortex 
and its exterior [as also suggested by the GNG vortex so- 
lution, cf. Eq. 1)23(1 ]. We also find that shock viscosity helps 
to reduce this problem. 

For the hot corona, we plot in Fig. gas density, gas 
velocity, dust density and dust velocity in the mid-plane 
after one orbit. The gas configuration is indistinguishable 
between different values of J?r s , indicating that there is 
very little back-reaction on the gas, so gas density and 
velocity field are only shown for fir s = 1.0. 

The vortex is evidently still in place, although the 
outer parts have been sheared away to form long tails in 
the direction of the shear. In Fig. 0gas density and dust 
density in the x-z plane are plotted. Again we see that the 
vortex is still in place after one orbit. We conclude from 
this that the hot corona vortex solution is valid even in 
the presence of shear. 

The status of the cylindrical vortex after one orbit is 
shown in Fig. |H1 In the mid-plane it seems to be more 
disrupted than the hot corona vortex, possibly due to the 
smaller density ratio between the vortex and its surround- 
ings. A vertical cut (not shown) reveals that the original 
stratification of the vortex is almost intact, which implies 
that the cylindrical vortex model is also a valid one. 

To test the lifetimes of vortices beyond one orbit we 
have run low-resolution simulations (64 x 128 x 64) for 
different values of the viscosity v. Since the major effect 
of viscosity takes place at the vortex boundary where the 
velocity (and thus the mass flux) is highest, we plot the 
maximum mass flux (pw) ma x in the box as a function of 
time measured in orbits in Fig. [5J The lifetime is obvi- 
ously very dependent on the viscosity, with lower viscosi- 
ties leading to longer-living vortices. A low viscosity, on 
the other hand, also leads to more chaotic development 
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Figure 8. Gas and dust in the mid-plane for the cylindri- 
cal vortex after one orbit, Qt s = 1.0. The vortex is close 
to breaking up: two arcs of matter are being expelled at 
the NE and SW corners. Note that the density contrast 
between the vortex and its surroundings is much smaller 
than for the hot corona vortex, and that the gas density 
colour scale is different from the hot corona's. 



in maximum mass flux, probably due to too high a mesh 
Reynolds number. 

7.4. Evolution of dust density 

Also shown in Figs and is the evolution of dust density. 
As expected, the most efficient dust-trapping occurs when 
dust and gas are semi-coupled with i?r s = 1.0, where the 
increase in dust density inside the vortex is from an initial 
In Pd/po = —4.18 to a maximum of In Pd/po = —1-34 after 
one orbit, a 17 times increase in density. A depletion of 
dust in the outer parts of the vortex is also seen. We be- 
lieve this is where the trapped dust has been taken from. 
There is a strong convergence in the dust velocity field 
both in the mid-plane and in the x-z plane. 

For Qt s = 0.1 dust is accelerated to match the gas ve- 
locity very quickly (compared to an orbit) and only a very 
modest increase in dust density is seen inside the vortex. 
This could be the slow settling of dust that was discussed 
in Sect. which should here happen on a time scale of 
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Figure 9. Maximum mass flux in the box as a function 
of time (measured in orbits) for different values of the 
viscosity v. A low viscosity obviously increases the lifetime 
of a vortex. On the other hand a low viscosity also leads 
to more chaotic development in (pw) ma x, indicating that 
the mesh Reynolds number may become too high. 



about 10/ J?. The velocity field of Fig. [7|does however not 
show any convergence inside the vortex. But according 
to Eq. Q34fl the inwards drift could happen with a veloc- 
ity of only u( ter V(W2) w T s f25 2 w 0.02 for ttr s = 0.1 
and e = 0.4, a contribution in peril of drowning com- 
pletely in the erratic motions present after one orbit, al- 
though it could still be present underneath it. A nar- 
row region of dust depletion is apparent around the edge 
of the vortex. This may be caused by the slow inwards 
drift of dust (its radius of ~ 0.1 is comparable to the ex- 
pected 0.02 • 2ir = 0.13), although it is not clear why the 
fir s = 10.0 vortex should have a region of similar deple- 
tion. In Fig. w e show gas and dust configurations of 
the J?t s = 0.1 vortex after 1.5 orbits. The dust-depleted 
region on the vortex boundary has obviously not become 
wider, nor has the dust density inside the vortex increased. 
Gas density has changed a lot. The vortex now has a pro- 
nounced tilt in the NW-SE direction, and the shear tails 
have grown more massive in density. The reason why we 
do not see any dust-trapping may then very well be that 
vortex dynamics completely dominates over the slow in- 
wards drift of dust, rendering the mechanism for trapping 
short stopping time dust particles useless. 

For a stopping time parameter of i?r s = 10.0 dust is 
so unaffected by gas that only a slight density increase is 
seen. The increase has occurred in a narrow region that 
extends from the vortex along the shear. This may be the 
dust-trapping mechanism that here takes place so slowly 
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Figure 10. Gas and dust density and velocity fields after 
1.5 orbits for Qt s = 0.1. Comparing with Fig. (plot 
number two from left) there is no apparent widening of 
the dust-depleted region on the vortex boundary, nor has 
the dust density increased inside the vortex. The colour 
scale is the same is in Fig. 0J 



that dust is sheared away faster than the vortex can trap 
it. 



7.5. Back-reaction on gas 

One very important issue regarding vortices is how they 
are affected by the dust they trap. Vertical dust-settling 
at the same time increases dust density in the mid-plane, 
and at one point dust density in the vortex is expected to 
become comparable to gas density. Then the back-reaction 
on gas becomes important in the dynamics of the gas, 
Eq. (JI3J. 

As mentioned in Sect. 17.31 the reason why gas develops 
similarly in our simulations for all values of stopping time 
is that the back-reaction on the gas from dust is negligible. 
This is caused by two things: a dust-to-gas ratio of 0.01 is 
not very high, and furthermore this ratio only applies in 
the disc, whereas we know that gas density in the vortex 
is very much higher, giving an even lower dust-to-gas ratio 
there. The back-reaction would eventually become much 
bigger if we allowed for vertical settling of dust. 








Figure 11. Gas configuration after half an orbit (left) and 
after a whole orbit (right) for a dust-to-gas ratio of the 
order of 1 inside the vortex (100 in the disc). Velocities are 
reduced by around a factor of 5 compared to the initial 
condition. Gas is forced to follow the converging dust, and 
that leads to the destruction of the vortex. The colour 
scale is the same is in Fig. 0] 



To test when the back-reaction from dust becomes im- 
portant we have run the J7r s = 1.0 and Qt s = 10.0 vortices 
with higher values of the dust-to-gas ratio. In Fig. 1111 we 
plot the gas configuration of the S7t s = 1.0 vortex when 
the dust-to-gas ratio is 100. At this point dust density 
becomes comparable to gas density inside the vortex. We 
see that the vortex is not only affected by dust drag, but 
that it is even almost completely destroyed after one or- 
bit. When dust converges inside the vortex, it will drag 
gas along with it, in this way destroying the vortex. 
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Figure 12. The maximum mass flux in the box after one 
orbit when varying the dust-to-gas ratio from 0.01 to 100, 
the latter corresponding to a dust-to-gas ratio of unity 
inside the vortex. Points are shown for both Qt s = 1.0 
and Qt s = 100.0. The best dust-trapper, Qt s — 1, is also 
the best at destroying the vortex, due to the trapped dust 
dragging gas along with it. 



The maximum mass flux in the box after one orbit 
as a function of dust-to-gas ratio in the disc is shown in 
Fig. ^] Already when the dust-to-gas ratio in the disc is 
10 (corresponding to around 0.1 in the vortex), the back- 
reaction becomes important. The f2r s = 1.0 vortex is best 
at destroying the vortex at all dust-to-gas ratios. 

8. Conclusions 

In this paper we have explored the suggestion of Barge & 
Sommeria ( 1995) of trapping dust particles in anticyclonic 
vortices. In recent years the vortex theory has gained in- 
creasing popularity as a site for planetesimal formation. 

To our knowledge this is the first time the dust- 
trapping mechanism has been explored in three dimen- 
sions and with a freely developing gas velocity field and 
density. We use two ways to model the vortices in three 
dimensions, and we have shown that both models sur- 
vive well for at least one orbit. It is conceivable that long 
life times are possible at yet higher resolution when the 
viscosity can be smaller. We have also shown the dust- 
trapping capabilities of vortices. For both vortex models 
dust-trapping is very efficient when dust and gas are semi- 
coupled, whereas too little or too much coupling gives 
no significant dust-trapping. The most efficiently trapped 
solid particles have a radius of 100 centimetres. This may 
be in some conflict with meteoritic evidence, where chon- 



drules (the building blocks of the most pristine meteorites) 
are typically up to 1 millimetre in radius. But meteorites 
on Earth are found to originate in the asteroid belt and 
should therefore not necessarily say anything about the 
make up of pristine material in the outer Solar System. 

Our two three-dimensional models both have a large 
density contrast between the vortex and the surrounding 
disc. A model with a realistic disc density profile could be 
obtained by using coordinates more natural to the vortex 
flow, such as elliptical coordinates. This would also require 
a version of the shearing sheet approximation adapted 
to new coordinates (for spherical coordinates there exists 
the shearing disc approximation of Klahr & Bodcnhcimer 
2003). Alternatively, and probably better, a global solu- 
tion to the Euler and continuity equations with a gradual 
transition from vortex to surrounding Kepler flow could be 
used. Such global solutions have been found to the Euler 
equation alone (Chavanis 1200.01 de la Fuente Marcos &; 
Barge 2 0010 , but none to our knowledge satisfy the conti- 
nuity equation. 

The role of vortices as high pressure regions seems 
hitherto unexplored. Short stopping time solid particles 
forced to match gas speed around the vortex will feel an 
additional pull inwards. The drift is similar to the vertical 
settling of dust towards the mid-plane. Although this hap- 
pens on quite a long time scale compared to the conven- 
tional dust-trapping time, it seems a viable mechanism for 
catching short stopping time dust particles, provided vor- 
tices live long enough. Unfortunately our simulations do 
not show this mechanism at work, but rather that vortex 
and gas dynamics completely dominates over this effect, 
even if it is present. 

It will be necessary to address the questions as to what 
causes such vortices in the first place. It seems plausible 
that anticyclonic vortices form as a relic after the disc tur- 
bulence in the disc has died out (Bracco et al. 119981 IT999). 
This idea is particularly attractive, because the absence of 
turbulence may be an important prerequisite for allowing 
dust to settle to the mid-plane. Other suggestions include 
baroclinic instabilities (Klahr & Bodenheimer, 2003), the 
anisotropic kinetic alpha-effect (von Rekowski et al.©99), 
and perhaps even the possibility of brown dwarf encounter 
(Willerding 2002 ). The original work by GNG was inspired 
by numerical results of the Papaloizou-Pringle instability, 
but this instability is not available in a thin Kepler disc 
with q = 3/2. Stability analysis in GNG suggested that 
the vortices suffer from numerous linear instabilities, so 
the survival of the vortices in 3D for at least one orbit 
seems surprising. In any case, it will be important to re- 
lax the need for implementing vortices as initial conditions 
and rather try to get them automatically under realistic 
conditions. 

The photometric observations of the protostar KH15D 
has been interpreted by Barge & Viton (2003) as a giant 
dusty vortex rotating at a distance of 0.2 AU from the 
star and covering 120° of the orbit. Such a large vortex 
is beyond the validity of the shearing sheet approxima- 
tion (where the curvilinearity of the disc is neglected) and 
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requires a global hydrodynamical disc simulation to exam- 
ine its stability and creation. If the vortex interpretation 
is correct, it could be the opening of a new era of observa- 
tional vortex research. The next generation of telescopes 
such as ALMA (see Wolf & Klahr !2UU2|l will be able to 
probe protoplanetary discs all over the sky for evidence of 
vortex activity. This will be the ultimate test of whether 
the research done so far in the field has been just an intel- 
lectual exercise, or if planets do indeed form partially as a 
result of long-lived vortices. This would make the vortex 
dust-trapping theory an important step in planet forma- 
tion. Perhaps we owe the existence of our own planet to 
vortices that were present in the solar nebular 4.6 billion 
years ago. 
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Appendix A: Derivation of vortex solution 

In this Appendix we show how to construct an elliptical 
velocity field and an appropriate corresponding enthalpy 
that satisfies Eqs. (|18() and (|19() . We start out by proposing 
that the velocity field Eqs. I|20[l - H22[) i s indeed a solution, 
provided that the aspect ratio and the angular velocity 
fulfil certain criteria. This velocity field is divergence-free 
and has no z-component. It can then be written as the 
curl of a stream function with only a z-component, 
u = Vxf. 

Given the velocity field u one can attempt to construct 
an enthalpy such that the flow is steady with du/dt = 0. 
This implies that 



Vh = f2 2 ('3x - z) - 2f2 x it - (u ■ V) 



(AT) 



The tidal force term and the vertical gravity term on the 
right hand side are obviously gradient terms. The same is 
true for the Coriolis term since 

/0\ / W/dy \ (W/dx\ 
kxu = ( I x I -dV/dx I = ( m/dy \ = W.(A.2) 

We are now left with (ignoring the integration constants 
for now) 

Vh = V(|r2 2 x 2 - \Q 2 z 2 - 2m) -(u-V)-u. (A.3) 

The advective term (u ■ V)u is a bit more tricky to put 
in gradient form. It can be done by applying the vector 
identity 

(u ■ V)u = V{\u 2 ) -u x (V x u) . (A.4) 

The first term on the right hand side comes out on gradient 
form. The second can be calculated by noting that 



V x u = oj — 






-d 2 1>/dx 2 - d 2 V/dy 2 



(A.5) 



where the vorticity lj of the flow is introduced. The sec- 
ond term on the right hand side of Eq. I|A.4(I can now be 
rewritten as 

dV/dy \ ( \ 
a (V \ u) = j -dtf/dx I x I I = -wVW . (A.6) 
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If the vorticity is constant and independent of the spatial 
coordinates, this can be written 



U X (V X It) = V(-LUlP) 



(A.7) 



The gradient of the enthalpy can now be written entirely 
as a sum of gradient terms, 



2„2 



. (A. 



This is easily integrated to give 



= m 2 x 2 



\Q 2 z 2 - (2f2 + \u 2 + const, (A.9) 



where all the integration constants are collected in just 
one constant. 

In order to calculate the enthalpy of the flow given by 
Eqs. I|20() - (|22|l we need to know also the stream function 
and the vorticity of the flow. The magnitude of the stream 
function is 



1 



and the vorticity has the magnitude 
d 2 <P 



dx 2 dy 



d 2 V /l 

= - - + e 



n' . 



(A.10) 



(A.ll) 



This is a constant, and Eq. (|A.9|I can therefore be ap- 
plied. To simplify the calculations we write the enthalpy 
ash(x,y,z) — h x (x) + h y (y) + h z (z) (as there are no mixed 
terms). For h x we get 



h x = {\fl 2 - Q-fl' + \fl' 2 )x 2 ee Bx 2 : 



similarly for h y , 

h y = {-Qef2' + \n l2 )y 2 

and for ft z , 

> 2 ~ 2 - Cz 2 . 



Ay 2 , 



(A.12) 



(A.13) 



(A.14) 



Now the enthalpy can be written as ft. = Bx 2 + Ay 2 + Cz 2 . 
This can be done for any angular velocity fl' of the vortex. 

An equilibrium flow solution must also have dp/dt = 
everywhere. The continuity equation Eq. I|19l) can be 
rewritten as 



dp 
di 



+ V • (pu) 



dp 

at 
dp 

dt 



•p(V-u) 
u ■ Vp = 



u ■ Vp 



0. 



(A.15) 



since the velocity field has V • u = 0. This means that the 
gradient of the density must everywhere be perpendicular 
to the flow, since u ■ Vp must be equal to zero. In the 
absence of an entropy gradient the gas is barotropic, so 
Vp and Vft are parallel, and therefore u ■ Vft = 0. This 
means that the contours of enthalpy must be ellipses with 
the same aspect ratio as the vortex. 



For the enthalpy given by h(x, y, z) — Bx 2 +Ay 2 
where the coefficients are specified in Eqs. (IA.12|I -(I 



"Oil 



and a velocity field given by Eqs. H2()|l - H22|l . this leads to 



A 
B 



-fief}' 



(A.16) 



which implies that 

ln 2 e 2 - n<tn' + h 



n' = 



V3ef2 



2 il 



ee af2. 



(A.17) 



This solution requires that ^ e < 1. Negative e in the 
same range in principle also give solutions, but the en- 
thalpy and velocity field do not change since e only enters 
as e 2 when fi 1 is inserted. The parameters B, A and C are 



B = -l[2 2 -- 



A 

C 



\n 2 



±n 2 



2V3 

7T = ? 

2V3 



(A.18) 

(A.19) 
(A.20) 



For 0.5 < e < 1 the coefficients B and A are positive, 
resulting in a region of low pressure with a clockwise ro- 
tation around it, contrary to the counter-clockwise rota- 
tion around low-pressures on Earth. As in the GNG paper 
we focus only on ^ e ^ 0.5. Here B and A are neg- 
ative, and the result is a high-pressure region. Defining 
5 2 = — + -7^=5 and requiring that the enthalpy van- 
ish on the vortex boundary gives 



h = U 2 n 2 {b 2 



\n 2 z 2 . 



(A.21) 



